Introduction à sf

• Site web de sf: Simple Features for R

• sf pour Simple Features

• Sortie en octobre 2016

• A pour but de rassembler les fonctionnalités d’anciens packages (sp, rgeos and rgdal) en un seul

• Facilite la manipulation de données spatiales, avec des objets simples.

• Tidy data: compatible avec la syntaxe pipe %>% et les opérateurs du tidyverse.

• Principal auteur et mainteneur : Edzer Pebesma (également auteur du package sp)

la structure des objets sf :

Importer / exporter des données

library(sf)

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

mtq <- read_sf("data/mtq/martinique.shp")
mtq <- st_read("data/mtq/martinique.shp")

## Reading layer `martinique' from data source `/home/comeetie/Projets/quantilille/lecture/data/mtq/martinique.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 23 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
## Projected CRS: WGS 84 / UTM zone 20N

write_sf(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)
st_write(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)

## Deleting layer `martinique' using driver `GPKG'
## Writing layer `martinique' to data source `data/mtq/martinique.gpkg' using driver `GPKG'
## Writing 34 features with 23 fields and geometry type Polygon.

Système de coordonées

Les projections/système de coordonées sont répertoriées grâce à un code le code epsg :

• lat/long : 4326 https://epsg.io/4326
• Lambert 93 : 2154 https://epsg.io/2154
• Pseudo-Mercator : 3857 https://epsg.io/3857
• Lambert azimuthal equal area : 3035 https://epsg.io/3035

Projection

Obtenir la projection en utilisant st_crs() (code epsg) et la modifier en utilisant st_transform().

st_crs(mtq)

## Coordinate Reference System:
##   User input: WGS 84 / UTM zone 20N 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 20N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 20N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-63,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     ID["EPSG",32620]]

mtq_4326 <- mtq %>% st_transform(4326)

Afficher les données

Affichage par défaut :

plot(mtq)

## Warning: plotting the first 10 out of 23 attributes; use max.plot = 23 to plot
## all

En ne gardant que la géométrie :

plot(st_geometry(mtq))

Extraire les centroïdes

mtq_c <- st_centroid(mtq)

## Warning in st_centroid.sf(mtq): st_centroid assumes attributes are constant over
## geometries of x

plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)

Matrice de distance

mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]

## Units: [m]
##           [,1]     [,2]      [,3]      [,4]      [,5]
## [1,]     0.000 35297.56  3091.501 12131.617 17136.310
## [2,] 35297.557     0.00 38332.602 25518.913 18605.249
## [3,]  3091.501 38332.60     0.000 15094.702 20226.198
## [4,] 12131.617 25518.91 15094.702     0.000  7177.011
## [5,] 17136.310 18605.25 20226.198  7177.011     0.000

Agrégation de polygones

Union simple :

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")

A partir d’une variable de regroupement :

library(dplyr)

## 
## Attachement du package : 'dplyr'

## Les objets suivants sont masqués depuis 'package:stats':
## 
##     filter, lag

## Les objets suivants sont masqués depuis 'package:base':
## 
##     intersect, setdiff, setequal, union

mtq_u2 <- mtq %>% 
  group_by(STATUT) %>% 
  summarize(P13_POP = sum(P13_POP))
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u2), add=T, lwd=2, border = "red", col=NA)

Zone tampon

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")

Intersection de polygones

# create a polygon
m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586), 
           c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)

st_intersection() extrait la partie de mtq qui s’intersecte avec le polygone créé.

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)

Compter les points dans des polygones

st_sample() crée des points aléatoires sur la carte de Paris.

pts <- st_sample(x = mtq, size = 50)
plot(st_geometry(mtq))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

st_interects() crée une liste de points dans chaque polygone.

inter <- st_intersects(mtq, pts)

Il ne reste plus qu’a compter.

mtq$nbpts <- sapply(X = inter, FUN = length)
plot(st_geometry(mtq))
# display munucipalities that intersect at least 2 point
plot(st_geometry(mtq[mtq$nbpts>2,]), col = "grey", add=TRUE)
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

Autre solution, faire une jointure spatiale et aggréger !

mtq_counts <- mtq %>% st_join(st_as_sf(pts)) %>% count(INSEE_COM)
plot(mtq_counts %>% select(n))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

Voronoi Polygons

google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)

mtq_v <- st_collection_extract(st_voronoi(x = st_union(mtq_c)))
mtq_v <- st_intersection(mtq_v, st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c)
plot(st_geometry(mtq_v), col='lightblue')

Autre traitements

• st_area(x)
• st_length(x)
• st_disjoint(x, y, sparse = FALSE)
• st_touches(x, y, sparse = FALSE)
• st_crosses(s, s, sparse = FALSE)
• st_within(x, y, sparse = FALSE)
• st_contains(x, y, sparse = FALSE)
• st_overlaps(x, y, sparse = FALSE)
• st_equals(x, y, sparse = FALSE)
• st_covers(x, y, sparse = FALSE)
• st_covered_by(x, y, sparse = FALSE)
• st_covered_by(y, y, sparse = FALSE)
• st_equals_exact(x, y,0.001, sparse = FALSE)
• …

Conversion

• st_cast
• st_collection_extract
• st_sf
• st_as_sf
• st_as_sfc

Autres packages

CRAN task views permet d’avoir des informations sur les packages du CRAN pertinents pour des tâches reliées à certains sujets.

CRAN Task View: Analysis of Spatial Data:

• Classes for spatial data
  
• Handling spatial data
  
• Reading and writing spatial data
  
• Visualisation
  
• Point pattern analysis
  
• Geostatistics
  
• Disease mapping and areal data analysis
  
• Spatial regression
  
• Ecological analysis

Utiliser osmdata

osmdata permet d’extraire des éléments de la base de données gratuite et open-source OpenStreetMap. Cela peut nous servir par exemple pour récupérer des élements d’habillage fleuve / routes ou d’autres informations. La requête suit la nomenclature osm sur base de clés / valeures. Vous pouvez utiliser tagingo pour explorer l’ensemble des clés et valeures utilisées par la communauté OSM.

library(osmdata)
# Récupérer les routes principales grâce à osm
bb      <- iris.75 %>% st_transform(4326) %>% st_bbox()
q       <- opq(bbox = bb,timeout = 180)
qm      <- add_osm_feature (q, key = 'highway',value = 'motorway', value_exact = FALSE)
qt      <- add_osm_feature (q, key = 'highway',value = 'trunk', value_exact = FALSE)
qp      <- add_osm_feature (q, key = 'highway',value = 'primary', value_exact = FALSE)

motorway<- osmdata_sf(qm)
trunk   <- osmdata_sf(qt)
primary <- osmdata_sf(qp)

roads    <- c(primary,trunk,motorway)$osm_lines %>% st_transform(st_crs(iris.75))
roads.geom = st_geometry(roads)

# Récupérer le shape de la seine 
qr <- q %>% 
  add_osm_feature (key = 'waterway') %>% 
  add_osm_feature(key = "name:fr", value = "La Seine")
river <- osmdata_sf(qr)

river.geom <- c(st_geometry(river$osm_lines),st_geometry(river$osm_multilines)) %>% st_transform(st_crs(iris.75))

# Export road and river layers to shapefile
st_write(roads.geom, dsn = "data/osmdata/road.shp")
st_write(river.geom, dsn = "data/osmdata/river.shp")

Utilisons ces données pour habiller un peu notre carte:

# bbox est utilisé pour centrer sur Paris
bb <- st_bbox(iris.75)
par(mar=c(0.2,0.2,1.4,0.2), bg="ivory")
plot(st_geometry(iris.75), col = "ivory", border="ivory3", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
plot(st_geometry(river.geom),col="#87cdde",lwd = 3,add=TRUE)
plot(accidents.2019.paris %>% filter(grav==3) , pch = 20, col = "orange", add=TRUE, cex = 1)

## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 3), pch = 20, :
## ignoring all but the first attribute

plot(accidents.2019.paris %>% filter(grav==2) , pch = 20, col = "red", add=TRUE, cex = 1)

## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 2), pch = 20, :
## ignoring all but the first attribute

Géocodage

Géocoder c’est passé d’une adresse à une position géographique. En france, la Base d’adresse national permet de faire ce travail efficacement, en R le package banR permet d’intéroger l’API de la BAN simplement ce package doit être installé via github et permet ensuite de géocoder une colonne d’adresse en batch en spécifiant la colone contenant les adresse et éventuellement une colone contenant un code insee pour faciliter et préciser la requête. Pour des adresses internationales il est possible d’utiliser tidygeocoder qui peut intéroger différente API gratuite ou payante et qui fonctionne de manière assez similaire.

# avec banR
# remote
library(banR)
geo_banR = accidents.2019.paris %>% 
  filter(catv %in% c("VAE","EDP à moteur")) %>% slice(1:10) %>% 
  geocode_tbl(adresse = voie,code_insee = com) %>% select(latitude,longitude) %>%
  st_as_sf(coords = c("longitude", "latitude"),crs = 4326, agr = "constant") %>%
  st_transform(2154)

## Writing tempfile to.../tmp/RtmpIrPGvn/file239a2d1957f8.csv

## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## Please use the `file` argument instead.

## If file is larger than 8 MB, it must be splitted
## Size is : 666 bytes

## SuccessOKSuccess: (200) OK

## New names:
## * geometry -> geometry...10
## * geometry -> geometry...13

library(tidygeocoder)
geo_tidygeocoder = accidents.2019.paris %>% 
  filter(catv %in% c("VAE","EDP à moteur")) %>% slice(1:10) %>%
  mutate(addr = paste(voie,", Paris, France")) %>%  
  geocode(addr,method="osm") %>% select(lat,long) %>%
  st_as_sf(coords = c("long", "lat"),crs = 4326, agr = "constant") %>%
  st_transform(2154)

st_distance(geo_banR,geo_tidygeocoder,by_element = TRUE)

## Units: [m]
##  [1] 950.52141 382.66079 122.00239 181.62935  47.91376 642.02833  41.90755
##  [8] 305.16959 664.76010 664.76010

Petite introduction de sémiologie graphique

Cartes interactives mapview

Les cartes interactives ne sont pas forcément très pertinentes pour représenter des informations géostatistiques.

En revanche, elles sont utiles pour explorer les bases de données. Voyons un exemple avec mapview concernant les accidents mortels à Paris en 2019.

#remotes::install_github("r-spatial/mapview")
library(mapview)
mapviewOptions(fgb=FALSE)
# construire une carte avec certaines options pour les cercles
# avec mapview la taille des cercles reste constante quel que soit le zoom. 
# grav = 2 : individus tués
mapview(accidents.2019.paris %>%
          filter(grav==2))

On customise un peu…

mapview(accidents.2019.paris %>%
          filter(grav==2),
        map.types = "Stamen.TonerLite", legend=FALSE,
        cex=5, col.regions="#217844", lwd=0, alpha=0.9)

On customise encore un peu plus…

mapview(accidents.2019.paris %>%
          filter(grav==2) %>%
          mutate(age=2019-an_nais),
        map.types = "Stamen.TonerLite", legend=FALSE,
        cex="age", zcol="sexe", lwd=0, alpha=0.9)

ggplot2

library(ggplot2)
ggplot(accidents.2019.paris)+geom_bar(aes(x=catv,group=sexe,fill=sexe))

catv_ol = accidents.2019.paris %>% st_drop_geometry %>% count(catv) %>% arrange(n) %>% pull(catv)
gg = accidents.2019.paris %>% mutate(catv_o = factor(catv,levels=catv_ol)) %>% filter(catv_o %in% tail(catv_ol,10))

ggplot()+geom_bar(data = gg,aes(y=catv_o,group=sexe,fill=sexe))+
  scale_fill_brewer("Sexe",palette="Set1")+
  theme_bw()+
  labs(title="Nombre d'accidentés par type de véhicule et sexe",
       subtitle="à Paris en 2019, pour les hommes et les femmes ",caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")

gg = accidents.2019.paris %>% 
  st_drop_geometry %>% 
  filter(catv %in% tail(catv_ol,9)) %>% 
  count(catv,lum,sexe) %>% 
  add_count(catv,wt=n,name="tot") %>% 
  mutate(prop = n/tot)

ggplot(gg)+geom_point(aes(y=lum,x=sexe,color=prop,size=prop))+
  facet_wrap(~catv)+
  scale_color_distiller(palette="Reds",direction=1)+
  labs(title="Part d'accidentés par type de véhicule et sexe",
       subtitle="à Paris en 2019. ",caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")

library(tidyr)
catv_sel = c("Bicyclette","VL seul","VAE","VU seul","EDP à moteur","Scooter < 50 cm3")

gg = accidents.2019.paris %>% 
  st_drop_geometry %>% 
  filter(catv %in% catv_sel) %>% 
  count(catv,sexe) %>% pivot_wider(names_from = "sexe",values_from="n")

ggplot(gg) + geom_segment(aes(x="Homme",y=Masculin,xend="Femme",yend=Féminin,color=catv)) 

ggplot(gg)+
  geom_segment(aes(x="Homme",y=Masculin,xend="Femme",yend=Féminin,color=catv))+
  geom_text(data=gg %>% filter(catv!="Scooter < 50 cm3"),aes(x="Homme",y=Masculin,label=catv,color=catv),hjust="left")+
  geom_text(data=gg %>% filter(catv!="VU seul"),aes(x="Femme",y=Féminin,label=catv,color=catv),hjust="right")+
  scale_color_discrete(guide="none")+theme_bw()+
  labs(title = "Nombre d'accidentés suivant le sexe et le type de véhicule", 
       subtitle="à Paris en 2019", caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")

Intégrer des données spatiales avec geom_sf

ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = river.geom, colour = "azure",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  geom_sf(data = iris.75 %>%  st_centroid(),
          aes(size= nbacc), colour="#E84923CC", show.legend = 'point') +
  scale_size(name = "Nombre d'accidents",
             breaks = c(1,10,100,200),
             range = c(0,5)) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "ivory",color=NA),
        plot.background = element_rect(fill = "ivory",color=NA)) +
  labs(title = "Nombre d'accidents de la route à Paris par iris",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",x="",y="")

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

library(RColorBrewer)

acc = iris.75 %>% 
  st_join(accidents.2019.paris) %>% 
  group_by(INSEE_COM,do_union=TRUE) %>% 
  summarize(nb_acc=n(),
            nb_vl=sum(if_else(catv=="VL seul",1,0),na.rm = TRUE),
            nb_edp=sum(if_else(catv=="EDP à moteur",1,0),na.rm = TRUE),
            nb_velo=sum(if_else(catv=="Bicyclette",1,0),na.rm = TRUE))

## `summarise()` has grouped output by 'INSEE_COM'. You can override using the `.groups` argument.

acc

## Simple feature collection with 20 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
## Projected CRS: RGF93 / Lambert-93
## # A tibble: 20 x 7
## # Groups:   INSEE_COM [20]
##    INSEE_COM do_union nb_acc nb_vl nb_edp nb_velo                       geometry
##    <chr>     <lgl>     <int> <dbl>  <dbl>   <dbl>                  <POLYGON [m]>
##  1 75101     TRUE        285   117      8      27 ((651932 6861789, 651927.1 68…
##  2 75102     TRUE        173    63     10      19 ((650718.1 6863519, 650706.1 …
##  3 75103     TRUE        202    63     10      21 ((652440 6862594, 652433.4 68…
##  4 75104     TRUE        277    80     25      41 ((651966.6 6861729, 651911.6 …
##  5 75105     TRUE        234    75     14      28 ((651960.7 6859931, 651886.2 …
##  6 75106     TRUE        273   112     16      30 ((650807.6 6860419, 650795.2 …
##  7 75107     TRUE        444   202     13      37 ((650112.5 6861125, 650080 68…
##  8 75108     TRUE        721   333     17      32 ((650208.6 6862824, 650181.1 …
##  9 75109     TRUE        320   118     11      41 ((650787.8 6863551, 650780.1 …
## 10 75110     TRUE        442   162     11      57 ((652193.8 6863596, 652167.8 …
## 11 75111     TRUE        509   171     20      69 ((653355.9 6863227, 653352.2 …
## 12 75112     TRUE       1000   441     13      77 ((655285.2 6858619, 655221.2 …
## 13 75113     TRUE        798   350     15      62 ((651662.8 6858694, 651667.5 …
## 14 75114     TRUE        704   266     19      49 ((650148.5 6858070, 650116.1 …
## 15 75115     TRUE        858   364     14      78 ((645885.5 6859571, 645942.5 …
## 16 75116     TRUE       1429   766     16      67 ((645783.9 6859590, 645571.2 …
## 17 75117     TRUE        903   414     12      47 ((648839.5 6864557, 648790.3 …
## 18 75118     TRUE        718   318     17      46 ((651443.5 6864903, 651436.5 …
## 19 75119     TRUE        855   386     16      65 ((653702.1 6865070, 653697.5 …
## 20 75120     TRUE        826   358     15      55 ((654559 6863398, 654461.5 68…

bks <- round(10*quantile(100*acc$nb_velo/acc$nb_acc, na.rm=TRUE, probs=seq(0,1,0.2)))/10
bks

##   0%  20%  40%  60%  80% 100% 
##  4.4  6.6  7.7  9.8 12.1 14.8

acc = acc %>% mutate(txaccvelo= 100*nb_velo/nb_acc,txaccvelo_cat = cut(txaccvelo,bks)) 


ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = acc, aes(fill = txaccvelo_cat)) +
  geom_sf(data = river.geom, colour = "#87cdde",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  scale_fill_brewer(name = "Part (En %)",palette = "Reds", na.value = "grey80") +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "ivory",color=NA),
        plot.background = element_rect(fill = "ivory",color=NA)) +
  labs(title = "Part des Accidentés à vélos",
       subtitle = "par arrondissement à Paris en 2019",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")